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We investigate sandpile models where the updating of unstable columns is done according to a 
stochastic rule. We examine the effect of introducing nonlocal relaxation mechanisms. We find 
that the models self-organize into critical states that belong to three different universality classes. 
The models with local relaxation rules belong to a known universality class that is characterized by 
an avalanche exponent r m 1.55, whereas the models with nonlocal relaxation rules belong to new 
universality classes characterized by exponents r w 1.35 and r ~ 1.63. We discuss the values of the 
exponents in terms of scaling relations and a mapping of the sandpile models to interface models. 
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The concept of self-organized criticality (SOC), pro- 
posed in a seminal paper by Bak, Tang and Wiesenfeld 
jjj, has made considerable impact on a number of fields 
in the natural and social sciences. The paradigm of SOC 
is an idealized sandpile where grains added to a pile dis- 
sipate their potential energy through avalanches with no 
characteristic scale Early experimental studies of 

real sandpiles lead to clear disagreement with the nu- 
merical models: bounded distributions of avalanche sizes 
were observed instead of the expect power law depen- 
dence |ig-Q. On the other hand, recent rice-pile exper- 
iments found power-law distributions of avalanche sizes 
jL5| and tracer transit times |l(| . These results sparked a 
renewed interest in the study of sandpile models [fl6|-p3[ . 

A class of sandpile models with stochastic toppling 
rules, which is denoted "rice-pile" models, was found to 
display SOC in one dimension with a power-law distri- 
bution of avalanche sizes characterized by the critical ex- 
ponent t w 1.55 |l6| p0[ . Recently, Ref. pij proposed a 
mapping of the model in to the motion of a linear 
interface through a disordered medium p5[ ] . For this uni- 
versality class, to which we will refer to as the local linear 
interface (LLI) universality class, the mapping allows the 
determination of all the exponents characterizing the dy- 
namics of the pile [ ^4|j26| , and shows that r is clearly dif- 
ferent from the mean-field value r = 3/2 (see, e.g., p7j). 
Several other SOC models have been conjectured to be 
in the LLI class p3| , p^ , p8|p9| ] . Thus, the question arises 
of what mechanisms lead to the emergence of the LLI 
universality class for one-dimensional stochastic sandpile 
models with a preferred direction. Previously, sandpile 
models in higher dimensions have been classified accord- 
ing to their degree of 'directedness' P,|30|]. 

In this paper we undertake an investigation of the 



mechanisms responsible for the emergence of a given uni- 
versality class for one-dimensional sandpile models with 
stochastic toppling rules, and we discuss the reasons for 
the robustness of the LLI universality class. To this 
end, we study a class of models with stochastic toppling 
rules. We use as our basic model the one we proposed 
in Ref. jl7|] and investigate the robustness of the criti- 
cal behavior upon modification of the toppling process. 
Specifically, we study stochastic variants of the models 
proposed in Q . 

We find that the local models belong to the LLI univer- 
sality class and have r ~ 1.55, while the nonlocal models 
belong to new universality classes with r « 1.35 and 
t ~ 1.63. In order to understand the new exponent val- 
ues, we discuss scaling relations fulfilled by the exponents 
and the mapping to interface models. 

First, we define the class of one-dimensional models: 
The system consists of a plate of length L, with a wall 
at i = and an open boundary at i = L + 1. The profile 
of the pile evolves through two mechanisms: deposition 
and relaxation. Deposition is always done at i = 1, and 
one grain at a time. During relaxation we look at all 
active columns of the pile: A column i of the pile is con- 
sidered active if, in the previous time step, it (i) received 
a grain from column i — 1, for the local models, or from 
i — j, j = 1, 2, ..JV, for the nonlocal models, (ii) toppled 
a grain to column i + 1, or (iii) column i + 1 toppled one 
grain to its right neighbor. If a column i is active and the 
local slope Shi = h(i) — h(i + 1) > Si, then with proba- 
bility p(Shi) several grains will be toppled from column 
i. Here, we study the case 



p{Shi) = min {1, g{5h l - Si)} . 



(1) 
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where g < 1 is a parameter M]. The number of grains 
N to be toppled is determined by either the limited (I) 
or unlimited (u) rule M: 



N 



No, 

5hi - Si, 



(0 
(u) 



(2) 



The toppled grains are then redistributed either accord- 
ing to the local (L) or nonlocal (N) rule 



+ 1) = h(i + 1) + N, (L) 
h(i + j) = h(i + j) + l, j = l,...,N, (N) 



(3) 



Grains toppled to columns i > L leave the system. We 
designate the models by the two letters that describe the 
rules applied: the local limited and unlimited models are 
referred to as LI and Lu, while the nonlocal limited and 
unlimited models are referred to as Nl and Nu. The pa- 
rameter p(Shi) describes friction between the grains and 
accounts for the fact that a large range of "stable" slopes 
was observed in the rice-pile experiment. The parameter 
g accounts for the effect of gravity on the packing config- 
urations. For a discussion of the physical interpretation 
of the models see also Refs. Jl7| , ^8[ . 

We study the models in the slowly driven limit, i.e., we 
take the rate of deposition to be slow enough that any 
avalanche, which might be started by a deposited parti- 
cle, will have ended before a new particle is deposited. 
The simulations of the models show that each system 
quickly enters a steady state which is characterized by 
varying avalanche sizes and a complex structure in time 
p7| . The size s of an avalanche can be defined in a num- 
ber of ways. First, we follow the definition of Ref. [fL5| , 
and calculate the size of an avalanche as the total poten- 
tial energy dissipated in between deposited particles. In 
Fig. [j] we show the avalanche distributions for the vari- 
ous models. We find that the probability density function 
can be well described by the power law form 



P(s,L) = s-Tf(s/L v ), 



(4) 



where r and v are critical universal exponents. (Alter- 
natively, we have that P(s,L) — L~@ f(sj L v ), where 
/3 = vt.) We used plots of consecutive slopes for dif- 
ferent system sizes to obtain the estimates for r (see ta- 
ble ||). Then, these values were used to collapse the data 
according to Eq. (||) and extract v. 

The results from our numerical simulations show that 
the LI and Lu models belong to the LLI universality class 
with r k, 1.55 and v ~ 2.25, in agreement with the re- 
sults of Refs. |jl(|[l7). On the other hand, for the non- 
local models we obtain values of the exponents that sig- 
nal the existence of new universality classes. The nonlo- 
cal limited (Nl) model is characterized by the exponents 
t = 1.35±0.05 and v = 1.55±0.05. The combined change 
of the number of toppled particles together with a non- 
local relaxation leads to yet a new universality class: For 
the nonlocal unlimited (Lu) model, we obtain the expo- 
nent values r = 1.63 ± 0.02 and v = 2.75 ± 0.05. We 
note that simple power laws, as in Eq. ([|), in all cases 
provide us with nice data collapses. This result should 



be contrasted with the investigation of similar rules for 
the (nonstochastic) models in Ref. Q where the results 
could not be described by a simple power law but instead 
required a multifractal scaling form. 

To further test our conclusion regarding the univer- 
sality classes for the rice-pile models, we study different 
definitions of avalanche size: When using the total num- 
ber of topplings, we find the same values for r and v as 
quoted above, which is due to the fact that a toppling 
event on average dissipates a fixed amount of potential 
energy (for a bounded distribution of slopes). For the 
lifetime T, we find that the probability density function 
is well described by the scaling form 



D(T,L) =T~ v g(T/L 



(5) 



with the exponent values listed in table | (see also Fig. ||). 
These results reassure us that the LI and Lu models in- 
deed belong to the LLI universality class, whereas the Nl 
and Nu models belong to new universality classes. 

Next, we discuss scaling relations that are obeyed by 
the critical exponents. It is well established that the av- 
erage avalanche size scales as 



(s) ~ L\ 



(6) 



with the value of the exponent q depending on how the 
pile is driven. Here we have q = 1 for all models (see, 
e.g., H and |^3| for other q values). Combining Eqs. (^) 
and we obtain the exponent relation 



v(2-T)=q, 



(7) 



which is fulfilled by all models. We have used (s) ~ 
L v { 2 - T ) ( m combination with (s 2 )) to obtain alternative 
estimates of the critical exponents, and we obtain val- 
ues in complete accordance with those listed in table |. 
Similar relations to Eqs. @— (R) can be derived for the 
lifetimes regarding the averages (T), (T 2 ), and then used 
to extract the critical exponents. Notice, however, that 
for the Nu model, (T) is a constant since we have y > 2. 
From the scaling of the probability densities (ji|) and (||) , 
we obtain 



<r-l)=a(y-l), 



(8) 



in agreement with our results (s ~ T", with u> = vjo). 
The above relations imply that there are only two inde- 
pendent exponents. 

As noted above, it was shown in Ref. |24|] that the rice- 
pile model in fl6| can be mapped to a linear interface 
model described by the continuum equation 



dH _ D d 2 H 



ot 



d.i 



n(x,H) 



(9) 



Here, H, which is obtained as h(x,t) = H(x — l,t) — 
H(x,t), counts the number of topplings of a given col- 
umn and rj(x, H) is a quenched "noise" which is related 
to the stochastic toppling probability. The rice-pile dy- 
namics imposes a driving of the interface at x — 0. The 
mapping predicts that v = 1+X, where \ is the so-called 
roughness exponent characterizing, e.g., the scaling of 
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the width of the interface with system size L (see, e.g., 
p2|). Numerically one has \ ~ 1-25 |^5| in excellent 
agreement with the value v k, 2.25. In addition, the cut- 
off exponent for the avalanche lifetime is a = z, where 
z is the so-called dynamic exponent which describes the 
propagation of correlations. 

The results presented here show that the LI and Lu 
models belong to the LLI universality class. This result 
can be easily understood from the mapping to the linear 
interface model: The toppling of several particles does 
not change the fact that the growth of the interface is still 
local. Since the surface tension term is the relevant term 
(in the renormalization group sense) it follows that the 
coarse-grained behavior of the models is governed by the 
interface equation @ as for the models in Refs. [ ^qjr^| . 

On the other hand, for the Nl and Nu models, the non- 
local toppling rules generate nonlocal growth which affect 
the interface motion. It has previously been shown that 
nonlocal interactions will in general lead to the emer- 
gence of new universality classes |33| , [34f which is con- 
firmed by our observations for the Nl and Nu models. 
The critical exponents for the interface equations corre- 
sponding to the Nl and Nu models are not known. We 
can, nevertheless, qualitatively understand the change in 
the values of the exponent r as follows: In general, we 
would expect that a nonlocal toppling rule would lead 
to a decrease of the value of r because more columns in 
the pile are perturbed at every time step, thus creating 
larger avalanches. This is indeed what happens for the 
Nl model. However, for the Nu model we observe an 
increase in the value of r. To understand this result, it is 
useful to look at the average slope of the pile and at the 
distribution of local slopes for the different models (see 
Fig. ||). For the limited models, the average slope remains 
practically the same upon the change of the number of 
particles toppled: 9.3 for the LI model and 9.1 for the 
Nl model. This means that the change from one (lo- 
cal) to many (nonlocal) columns being "actived" works 
as we described above, i.e., it leads to an lower value of 
r. On the other hand, for the unlimited models, a large 
change in the average slope is observed when we change 
the number of particles toppled: 7.8 for the Lu model 
and 6.2 for the Nu model. This implies that, on average, 
only a few particles are toppled from unstable columns 
for the Nu model, but for the Lu model more particles 
are toppled and this leads to an higher likelihood of big 
avalanches. As a result, we conclude that the Nu model 
should have a higher value of r than the Lu model, as is 
indeed observed. 

In summary, we study a class of rice-pile models which 
in a simple way model some of the physical features of 
the experiments in Refs. We find that for local 

relaxation rules the models belong to the LLI universality 
class which is characterized by the exponents r k, 1.55 
and y « 1.87 (models LI and Lu). On the other hand, 
for models with nonlocal relaxation rules we obtain more 
complex dynamics resulting in new universality classes: 
For the Nl model we obtain r « 1.35 and y w 1.60, while 
for the Nu model we obtain r w 1.63 and y w 2.20. Our 
results show that nonlocal rules can increase or decrease 



the value of r, thus opening for the possibility that the 
rice-pile experiment in Ref. |15[| can be explained by some 
extension of nonlocal sandpile models by incorporating in 
detail the rice-grain dynamics. 
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Lifetime 

FIG. 2. Log-log plot of the probability density for the 
lifetime of avalanches. For greater clarity, the data for the 
Lu, Nl, and Nu was divided by factors of 10, 10 3 and 10 4 , 
respectively. The same values were used for the parameters 
as in Fig. hj. The slopes of the straight lines correspond to the 
exponent values from Table Nl 
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FIG. 1. Log- log plot of the probability density of the dis- 
sipated potential energy s during an avalanche. For greater 
clarity, the data for the Lu, Nl, and Nu was divided by factors 
of 10\ 10 2 and 10 3 , respectively. P(s) was multiplied by s to 
make visually more clear the change in the exponent r. For 
each model, we performed runs with the order of 10 7 grains 
deposited. The results shown are for L — 1600, g = 1/8, 
Si = 6, and No = 2 (for the limited models). Other values 
of Si (ranging from 1 to 6), of 1/g (ranging from 4 to 8), or 
of No (ranging from 2 to 4) were also investigated without 
any observable change in the estimate of the exponents. The 
slopes of the straight lines correspond to the exponent values 
from Table |. For the Nl model the bump before the cutoff 
leads to finite size corrections, and we carried out simulations 
on systems up to size L = 20000 in order to obtain more 
accurate estimates for the exponents. 
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FIG. 3. Linear-log plot of the probability of the local 
slopes in the steady state. The parameters of the runs are 
the same as in Fig. hi except that the system size is L = 400. 
It is visually apparent, that while the limited models lead to 
truncated Gaussian distributions, the unlimited models lead 
to more complicated forms. In particular, the Nu model leads 
to a nearly exponential decay of the probability of finding 
slopes larger than S\ (which is related to the increase in the 
value of r for this model). 



TABLE I. Critical exponents for the four models studied. 
The definition of the critical exponents and classification of 
the models are given in the text. The data strongly suggest 
that the local models, LI and Lu, belong to the LLI univer- 
sality class whereas the nonlocal models, Un and Ln, belong 
to new universality classes. 



Model 


r 


V 


V 


a 


LI 


1.55 ±0.02 


2.24 ±0.03 


1.83 ±0.04 


1.42 ± 0.03 


Lu 


1.56 ±0.02 


2.26 ±0.03 


1.91 ±0.04 


1.37 ±0.03 


Nl 


1.35 ±0.05 


1.55 ±0.05 


1.60 ±0.05 


0.95 ± 0.05 


Nu 


1.63 ±0.02 


2.75 ±0.05 


2.20 ± 0.04 


1.50 ±0.04 



4 



